use germany-pgd , replace

  /*******************************************************************************************/
  /* Regress positive/negative evaluations on religiosity, secularism, denomination, gender, */
  /* partisanship, medical background, experience with disability/childlessness,             */
  /* number of kids, age, education, east                                                    */
  /*******************************************************************************************/


  svy: reg pgdpositive c.religiosity c.sepp i.male i.denomination i.partisanship i.medicalbackground  i.exphandicap i.expchildless i.catkid  i.agecat i.educ i.ostsoz  
  est store pgdpositive

svy: reg pgdnegative c.religiosity c.sepp i.male i.denomination i.partisanship i.medicalbackground  i.exphandicap i.expchildless i.catkid  i.agecat i.educ i.ostsoz  
  est store pgdnegative


  /*******************************************************************************************/
  /* OGLM models                                                                              */
  /*                                                                                         */
  /*                                                                                         */
  /*******************************************************************************************/


  svy: oglm t2v1 c.pgdpositive##c.pgdnegative, het(c.pgdpositive##c.pgdnegative)
  est store fulloglm1

  * store coefficients for further analysis in R
  * parmest , saving(oglm1 , replace)

  svy: oglm t2v2 c.pgdpositive##c.pgdnegative , het(c.pgdpositive##c.pgdnegative)
  est store fulloglm2

  * store coefficients for further analysis in R
  * parmest , saving(oglm2 , replace)

  /****************************/
  /* store predictive margins */
  /****************************/

    * the idea is to store the margins and there CIs for three scenarios
  * these are then plotted with ggplot

  * these are the three modes -> scenarios
  * 1  1.125 -0.500
  * 2  0.750  0.375
  * 3 -0.375  0.375
  * we do this for each scenario X group combination

  * Group I
  est restore fulloglm1
  margins , at(pgdpositive=1.125 pgdnegative=-0.5) post
  *  parmest , saving(oglmmargins-1-1 , replace) idstr("I") idnum(1)

  est restore fulloglm1
  margins , at(pgdpositive=0.750 pgdnegative=0.375) post
  *  parmest , saving(oglmmargins-1-2 , replace) idstr("I") idnum(2)

  est restore fulloglm1
  margins , at(pgdpositive=-0.375 pgdnegative=0.375) post
  *  parmest , saving(oglmmargins-1-3 , replace) idstr("I") idnum(3)

  * Group II
  est restore fulloglm2
  margins , at(pgdpositive=1.125 pgdnegative=-0.5) post
  *  parmest , saving(oglmmargins-2-1 , replace) idstr("II") idnum(1)

  est restore fulloglm2
  margins , at(pgdpositive=0.750 pgdnegative=0.375) post
  *  parmest , saving(oglmmargins-2-2 , replace) idstr("II") idnum(2)

  est restore fulloglm2
  margins , at(pgdpositive=-0.375 pgdnegative=0.375) post
  * parmest , saving(oglmmargins-2-3 , replace) idstr("II") idnum(3)

* Additional analyses

  display "Positive attitude, group II, combined support for options 3-5"
  est restore fulloglm2
  margins, at(pgdpositive=1.125 pgdnegative=-0.5) expression(predict(pr outcome(3)) + predict(pr outcome(4)) + predict(pr outcome(5)))


  display "Positive attitude, group I, support for option 3"
  est restore fulloglm1
  margins, at(pgdpositive=1.125 pgdnegative=-0.5) expression(predict(pr outcome(3)) )

  display "Negative attitude, group I: show upper/lower limit for combination of 1+2"
  margins, at(pgdpositive=-0.375 pgdnegative=0.375) expression(predict(pr outcome(1)) + predict(pr outcome(2)) )
